home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byte0687.arc / POLYFIT.BAS < prev    next >
Encoding:
BASIC Source File  |  1987-03-05  |  4.7 KB  |  115 lines
  1. 1000 LN=1000: LD=11: REM LN=Max data points; LD=highest degree+1
  2. 1010 DEF FNMI(X,Y)=(X<Y)*(-X) + (Y<=X)*(-Y)
  3. 1020 N=0: M=0: S2=0: R2=0: MF=0
  4. 1030 S1=0: S2=0: S3=0: S4=0: P1=0: P2=0: P3=0: I=0: J=0: J1=0: K=0: VR=0: MM=0: WT=0: P=0
  5. 1040 DIM X(LN),Y(LN),W(LN),C(LD)
  6. 1050 DIM D1(LD),D2(LD),D3(LD),D4(LD),D5(LD),D6(LD)
  7. 1060 GOTO 1450
  8. 1070 IF MF>0 AND M>MM THEN J1=MM+1: MM=M: GOTO 1130
  9. 1080 J1=1: MM=M: S1=0: S2=0: S3=0: S4=0
  10. 1090 FOR I=1 TO N: WT=W(I)
  11. 1100 S1=S1+WT*X(I): S2=S2+WT: S3=S3+WT*Y(I): S4=S4+WT*Y(I)*Y(I)
  12. 1110 NEXT I
  13. 1120 D4(1)=S1/S2: D5(1)=0: D6(1)=S3/S2: D1(1)=0: D2(1)=1: VR=S4-S3*D6(1)
  14. 1130 FOR J=J1 TO MM: S1=0: S2=0: S3=0: S4=0
  15. 1140 FOR I=1 TO N: P1=0: P2=1
  16. 1150 FOR K=1 TO J: P=P2: P2=(X(I)-D4(K))*P2-D5(K)*P1: P1=P: NEXT K
  17. 1160 WT=W(I): P=WT*P2*P2
  18. 1170 S1=S1+P*X(I): S2=S2+P: S3=S3+WT*P1*P1: S4=S4+WT*Y(I)*P2: NEXT I
  19. 1180 D4(J+1)=S1/S2: D5(J+1)=S2/S3: D6(J+1)=S4/S2: D3(1)=-D4(J)*D2(1)-D5(J)*D1(1)
  20. 1190 IF J<4 THEN 1210
  21. 1200 FOR K=2 TO J-2: D3(K)=D2(K-1)-D4(J)*D2(K)-D5(J)*D1(K): NEXT K
  22. 1210 IF J>2 THEN D3(J-1)=D2(J-2)-D4(J)*D2(J-1)-D5(J)
  23. 1220 IF J>1 THEN D3(J)=D2(J-1)-D4(J)
  24. 1230 FOR K=1 TO J: D1(K)=D2(K): D2(K)=D3(K): D6(K)=D6(K)+D3(K)*D6(J+1): NEXT K
  25. 1240 NEXT J
  26. 1250 FOR J=1 TO M+1: C(J)=D6(M+2-J): NEXT J
  27. 1260 P2=0: FOR I=1 TO N: P=C(1)
  28. 1270 FOR J=1 TO M: P=P*X(I)+C(J+1): NEXT J
  29. 1280 P=P-Y(I): P2=P2+W(I)*P*P: NEXT I
  30. 1290 S2=0: IF N>M+1 THEN S2=P2/(N-M-1)
  31. 1300 R2=1: IF VR<>0 THEN R2=1-P2/VR: IF R2<0 THEN R2=0
  32. 1310 RETURN
  33. 1320 REM GOSUB 30 calls the subroutine
  34. 1330 REM Input: 
  35. 1340 REM N=# of data points
  36. 1350 REM X()=X-coordinates of the data points
  37. 1360 REM Y()=Y-coordinates of the data points
  38. 1370 REM W()=Weighting factors of the data points
  39. 1380 REM M=Degree of polynomial
  40. 1390 REM MF=0 if new data, MF=1 if old data but higher degree
  41. 1400 REM
  42. 1410 REM Output: 
  43. 1420 REM C=Array of M+1 coefficients
  44. 1430 REM S2=Residual variance
  45. 1440 REM R2=coefficient of determination
  46. 1450 CLS
  47. 1460 PRINT "Polyfit. Copyright (C) 1986 by William G. Hood"
  48. 1470 PRINT
  49. 1480 PRINT "This program finds the coefficients of the nth degree polynomial"
  50. 1490 PRINT: PRINT "y = c(1)*x^n + c(2)*x^(n-1) + ... + c(n)*x + c(n+1)
  51. 1500 PRINT: PRINT "that fits a set of data points in a least-squares sense."
  52. 1510 PRINT "Each data point must consist of an X value, a Y value, and an
  53. 1520 PRINT "optional weight, separated by commas and terminated by a return.
  54. 1530 PRINT "Data are read from a disk file until the eof is reached, or a
  55. 1540 PRINT "specified number of data points are read in from the keyboard."
  56. 1550 PRINT
  57. 1560 LINE INPUT "Name the data file (null line=keyboard): "; FI$
  58. 1570 PRINT "Is the data weighted (Y/N, null line=No)? ";
  59. 1580 W$="": INPUT W$
  60. 1590 IF W$<>"Y" AND W$<>"y" THEN W$="N"
  61. 1600 IF FI$<>NU$ THEN 2000
  62. 1610 PRINT "Keyboard data entry"
  63. 1620 PRINT "How many data points ( 2 <= n <=";LN;")";: INPUT N
  64. 1630 IF N<2 OR N>LN THEN 1620
  65. 1640 FOR K=1 TO N
  66. 1650 IF W$<>"N" THEN INPUT "x,y,w";X(K),Y(K),W(K): W(K)=ABS(W(K))
  67. 1660 IF W$="N" THEN INPUT "x,y"; X(K),Y(K): W(K)=1
  68. 1670 NEXT K
  69. 1680 PM=FNMI(LD-1,N-1)
  70. 1690 PRINT "Degree of polynomial ( 1<= m <= "; PM;")";: INPUT M: M=INT(M)
  71. 1700 IF M<1 OR M>PM THEN 1690
  72. 1710 GOSUB 1070
  73. 1720 PRINT: PRINT "Coefficients (constant term last): ": K=0
  74. 1730 FOR J=1 TO M+1: PRINT TAB(K) C(J);: K=K+20: IF K>20 THEN K=0: PRINT
  75. 1740 NEXT J
  76. 1750 PRINT: PRINT "Residual variance: ";S2
  77. 1760 R2=INT(1000000!*R2+.5)/1000000!
  78. 1770 PRINT: PRINT"Coefficient of determination: ";R2
  79. 1780 PRINT "Try another degree (Y/N, null line=No)";: A$="": INPUT A$
  80. 1790 IF A$="Y" OR A$="y" THEN 1680
  81. 1800 PRINT "Print a table of the data (Y/N, null line=No)";: A$="": INPUT A$
  82. 1810 IF A$<>"Y" AND A$<>"y" THEN 1870
  83. 1820 PRINT: PRINT"Degree of p(x): ";M
  84. 1830 PRINT: PRINT " X"; TAB(11) "Y"; TAB(25) "p(x)": PRINT
  85. 1840 FOR I=1 TO N: P=C(1)
  86. 1850 FOR K=1 TO M: P=P*X(I)+C(K+1): NEXT K
  87. 1860 PRINT X(I);TAB(10);Y(I);TAB(24);P: NEXT I
  88. 1870 IF F$<>NU$ THEN END
  89. 1880 PRINT: PRINT "Save the data (Y/N, null line=No)? ";: A$="": INPUT A$
  90. 1890 IF A$<>"Y" AND A$<>"y" THEN END
  91. 1900 PRINT: LINE INPUT "Name of output file ";FO$
  92. 1910 IF FO$=NU$ THEN END
  93. 1920 PRINT "Writing to ";FO$
  94. 1930 OPEN FO$ FOR OUTPUT AS 1
  95. 1940 FOR I=1 TO N
  96. 1950 IF W$<>"N" THEN PRINT#1, X(I);","; Y(I);","; W(I)
  97. 1960 IF W$="N" THEN PRINT#1, X(I);","; Y(I)
  98. 1970 NEXT I
  99. 1980 CLOSE 1
  100. 1990 END
  101. 2000 PRINT "Reading from ";FI$
  102. 2010 OPEN FI$ FOR INPUT AS 1
  103. 2020 N=0
  104. 2030 PRINT "  X" TAB(15) "Y"; TAB(28) "W"; : PRINT
  105. 2040 IF EOF(1) OR N=LN THEN 2100
  106. 2050 N=N+1
  107. 2060 IF W$<>"N" THEN INPUT#1, X(N), Y(N), W(N): W(N)=ABS(W(N))
  108. 2070 IF W$="N" THEN INPUT#1, X(N), Y(N): W(N)=1
  109. 2080 PRINT X(N) TAB(14) Y(N) TAB(28) W(N)
  110. 2090 GOTO 2040
  111. 2100 CLOSE 1
  112. 2110 PRINT: PRINT"File contained "; N;" data points."
  113. 2120 IF N<2 THEN PRINT "Too few data points.": END
  114. 2130 GOTO 1680
  115.